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ABSTRACT 

This paper presents detailed description of a novel CFD 
procedure and comparison of its solution results to that 
obtained by other available CFD codes as well as actual flight 
and wind tunnel test data pertaining to the GUI aircraft, 
currently undergoing flight testing at AFRC. 

INTRODUCTION 

Two in-house 1 software as well as a number of 
commercially 2 " 6 available CFD codes were used to analyze the 
problem, for comparison purposes. In this process both finite 
volume and finite element discretization were used for Euler 
and Navier-Stokes simulations. Both unstructured and 
structured grids were employed, as appropriate and solutions 
were derived for Mach 0.701 and angle of attack a = 3.92 
degree. 

Extensive flight tests were performed for validation purposes. 
Also these tests were complimented with detailed wind tunnel 
simulations. All such test results are compared with the 
numerical solution data obtained by the various CFD codes. 
Associated finite difference 7 and finite volume 8 - 9 techniques are 
well described in the literature 10,11 . The finite element 
technique 12 for the discretization of fluid flow employs 
unstructured mesh and is based on a Taylor-Galerkin 
procedure 13 ' 15 . 

A description of the finite element fluids solver is 
presented in some detail. It pertains to the solution of viscous 
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flow represented by the Navier-Stokes formulation. An 
unstructured grid is used for domain decomposition. 

The one equation model (Ref. 16) has been adapted for 
turbulence modelling. In this process both the viscous stresses 
pertaining to the linear viscous flow and the flux in the energy 
equation, are duly modified. 

It is then followed by detailed results of analyses which are 
next compared with actual flight test and wind tunnel 
simulation results. These results indicate that most CFD 
solutions compare reasonably well with the test data. The FE 
solutions in particular prove to be efficient and accurate and the 
related software are available for public use. 

Finally, some summarizations and discussions of the 
current effort is given in the ‘Concluding Remarks’ section. 


NOMENCLATURE 


AFRC 

- 

Armstrong Flight Research Center 

CFD 

- 

Computational Fluid Dynamics 

FE 

- 

Finite Element 

At 

- 

time step 

P 

- 

Density 
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- 

Dynamic viscosity 

a 

- 

Viscous stress tensor 

u 

- 

free stream velocity 

E 

- 

Total energy 

a 

- 

Shape function 

V 

= 

Conservation variable 

/ 

= 

Convection 

9 

= 

Diffusion 


1 



k 

P 

M 

K 

Re 

Pi- 


Thermal conductivity 
Pressure 
Mass matrix 
Convection matrix 
Reynolds number 
Prandtl number 


PROCEDURE 


The Navier-Stokes equation can be written as 

i =1,2,3 


— + — + — = 0 

dt dxi dxi 
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in which the conservation variables, flux, and body force 
column vectors, as well as the viscous stress are defined as 

v — [p puj pE] T , j = 1,2,3 (2) 
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where u L are velocity components in the x i coordinate system; 
p, p, and E are the density, pressure, and total energy 
respectively; p is the dynamic viscosity; k is the thermal 
conductivity, the heat flux qj being— kdT/dXj-, T is the 
temperature; f b represents the body forces. 

The preceding equations are nondimensionalised for 
numerical calculations. In this process the governing equations 
remain in the same form excepting g , , which becomes 

9j = [0 <J tJ ( UiCTij - qj)] T (7) 

and also the viscous stress tensor and heat flux take the 
following form: 


( 8 ) 


in which the Reynolds number is defined as Re — u m L/v m \ 
v oo — Moo/Poo is termed the kinematic viscosity; Pr is the 
Prandtl number, Pr—v m /a m , with = k/(p m c p ') is the 
thermal diffusivity. 

The Taylor’s expansion of the solution v(x, t) in the time 
domain, neglecting second order term and body forces, yields 
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in which Av — v(t + At) — v(t). Applying Galerkin’s spatial 
idealization v = av, v being the nodal values and a the shape 
functions vector, the flow equation can be expressed as 7 

MAv = -At + K} v - A t(J 1 + / 2 ) + AtR + 

At [K a + f a \ (10) 

in which M is the consistent mass matrix, K the convection 
matrix, / 1; f 2 the pressure matrices, K n the second-order 
matrix that includes viscous and heat flux effects, and f n the 
boundary integral matrix from second-order terms. Then, 
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In these equations, p ; , u h e, are the average values; e 1 — 
[0100 u 1 ] T , e 2 = [0010 u 2 ] t , e 3 = [0 0 0 1 u 3 ] T , R is the 
artificial dissipation, and m 1 — m 2 — m 3 = [0 0 0 0 1\ T . 
Turbulence terms are included by modifying the viscous 
effects. 

A novel two-step solution procedure 17 is adopted for the 
flow equation, the inviscid solution being augmented with the 
viscous term and stabilized with artificial dissipation terms. 
Assuming, 

Av = v n +i - v n (12) 

then. 
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Let 
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M + = D + 


+ M' + (17) 

the matrix D + having diagonal elements. Equation (8) may then 
be solved as follows. 

Step 1 : Form 


[D + ]v n+1 = [M_]v n - [M' + ]v n+1 + AtR 


(18) 


n+1 


Step 3: If || | 
Step 4: If Hv^l 


Step 2: Solve v n+1 iteratively 

- (i+1) = [D + ]<■-*> {[M_ ]v n - [M\ ]u® j + A t(R + R + 
Ka + fa )} (19) 

A EPSl||v® 1 || go to Step 2. 

A EPS2||v ) ( I I | 1 || go to Step 1. 

Step 5: Repeat Steps 1 to 4 NITER times until desired 
convergence is achieved, that is until v n+1 ~v n ; EPS1 and 
EPS2 are suitable convergence criteria factors, specified by the 
users. 

The iterative process in Step 2 requires a small number of 
steps, usually 1, and achieves a stable, convergent solution. 

In regions of high pressure gradients, artificial dissipation term 
is applied to prevent oscillations near discontinuities. This is 
implemented by incorporating pressure-switched diffusion 
coefficients as appropriate. Thus, 

C s S„ 


R = 


At 


-Ml 1 


[M c - M L ]v n 


( 20 ) 


in which C s is a shock capturing constant, S e is the averaged 
element value of the nodal pressure switch defined as 
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Z(|Pi -Pj\) 

and M c and M L are the consistent and lumped mass matrices 
respectively; l is the node under consideration and j are the 
nodes connected to i. 

To obtain the viscous components, ct l j in Eq. (4) is written 
as 


_ lg . n fdu t duj \ 

3 Re dxi 0 Re\dXj dxi) 

iffusioi 

9 i = (o a t 


( 22 ) 


and the diffusion flux of the Navier-Stokes equation being 
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i = 1,2,3; j = 1,2,3 (23) 

|i is the nondimensional viscosity term, whereas Rc and Pr are 
the Reynolds and Prandtl numbers, respectively. Next 
components of dgj dx L are evaluated term by term and then 
discretized by Galerkin approximation. 

This procedure is adopted in the STARS-CFDSOL code 1 
that enables effective solution of the Naviar-Stokes equation in 
most flight regimes. 


NUMERICAL AND TEST RESULTS 


right wing section was used for CFD analysis. Total number of 
CFD mesh using triangular element on wing surface is 31k for 
coarse mesh and 59k for finer mesh. Total number of 3-D CFD 
mesh using tetrahedron element in aerodynamic domain is 1 .2m 
for coarse mesh and 2.8m for finer mesh. 

. The flight condition was for Mach 0.7 and angle of attack a = 
3.92 degrees. Table 2 provides extensive description of relevant 
analyses hardware employed for each of the participative code 
and solution CPU time for a converged solution. The STARS 
has two solution option modules, namely CFDSOL and MG 
and both appear to be competitive in terms of solution time, 
accuracy, grid size and CPU numbers. 

Figures 4 to 6 depict pressure (C p ) distribution around the 
wing airfoil cross section at the wing 368.3 cm, 584.2 cm and 
1003.5 cm span wise locations. Further, the wind tunnel and 
actual flight test results are also shown for comparison and 
validation purpose. Due to the proprietary nature of the wind 
tunnel and fight test data, actual scales on the figures cannot be 
shown. Each of the codes shows reasonable correlation; 
solution of the CFDSOL and MG codes appear to be rather 
close to the two test results. 

Figure 7 depicts the C p distribution along the airfoil at 
different span locations. 


Accuracy 18 of the STARS CFD code was verified pertaining 
to the Hyper-X flight vehicle, carrying the X-43 vehicle for 
subsequent hypersonic flight at Mach 5.0 and 7.0. Table 1 
provides such a comparison of computational results and actual 
flight test data at various sensor locations; these data pertain to 
the ascent state of Hyper-X at Mach 0.9 and an altitude of 
22,500 ft. Figure 1 provides a graphical depiction of 
comparison of the two sets of results, signifying accuracy of the 
relevant procedures. Also Table 1 shows the numerical values 
of flight test and computed aerodynamic pressures; excellent 
correlation is observed for primary data values; the last three 
values in the Table are comparatively small and hence prone to 
measurement inaccuracy. This code was next used, along with a 
variety of existing commercially available programs, to solve a 
practical project problem. The results of which were also 
compared to that obtained by actual flight and wind tunnel 
tests. 

The Gulfstream GUI airplane (Gulfstream Aerospace 
Corporation, Savannah, Georgia), currently undergoing flight 
tests 19 at NASA AFRC, was chosen as the example problem for 
verification purposes. The GUI business jet as shown in Figure 
2 is being modified and instrumented by NASA's Armstrong 
Flight Research Center to serve as a test bed for a variety of 
flight research experiments, in support of the Environmentally 
Responsible Aviation (ERA) project. The twin-turbofan aircraft 
provides long-term capability for efficient testing of subsonic 
flight experiments for NASA, the U.S. Air Force, other 
government agencies, academia, and private industry. 

The wing span of the GUI aircraft is 23.7226 meter with 
sweep angle 27.66 degree. The airfoil section is a NACA 0012 
modification. The aerodynamic model of the GUI wing used in 
the CFDSOL and MG solutions is shown in Figure 3; only the 


CONCLUDING REMARKS 

The paper presents detailed comparison of solutions of the 
GUI aircraft wing obtained by a number of commercially 
available CFD codes as well as two AFRC in-house codes that 
use a finite element fluids discretization employing 
unstructured grids; related formulations of the novel CFDSOL 
code are also presented in detail. Importantly these solutions 
are compared with actual flight and also wind tunnel test data. 
Each of the codes shows reasonable correlation; solution of the 
CFDSOL and MG codes appear to be rather close to the two 
test results, particularly around the leading edge; further, use of 
a single CPU to derive solutions testifies to their cost 
effectiveness. 
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Table 1 Comparison of computed and flight test measured pressure 
data for the Hyper- X/X-43 vehicle 


Sensor 

point 


Pressure, Mpa 


Flight test 

CFD computed 

Percent 

difference 

001 

0.01165 

0.01193 

2.34 

003 

0.01227 

0.01164 

6.12 

007 

-0.00167 

-0.00096 

42.12 

085 

-0.00108 

-0.00268 

147.99 

090 

0.00048 

-0.00055 

2.56 


Table 2 CFD Solvers Comparison 


CFD Solver 

Flow Equation 

Platform 

No. of 
CPU 

Total CPU time 

Grid Size 

Note 

STARCCM+ 

RANS, finite 
volume, K- 
omega SST 
turbulence 

Cluster 

~80 

6hr, 40min (533 
cpu hours) - 
3000 iterations 

7.2M 

polyhedra/prismatic 
for half model 
without T-tail 

number of 
processers is 
an estimate, 
and the time 
is an estimate 
for that 
number of 
processors 

STARS (MG) 

Euler, finite 
element 

Dell M620 8GB 
Ram, 64 bit 

1 Intel 
Core i7 
@2.67 
GHz 

2.8 hr, (100 
steps, 25 inner 
cycles) 

1.2 M Tetrahedrons 
for wing only 


STARS 

(CFDSOL) 

Full N-S, finite 
element 

Dell M620 8GB 
Ram, 64 bit 

1 Intel 
Core i7 
@2.67 
GHz 

13.8 hr (10000 
steps) 

2.8 M Tetrahedrons 
for wing only 


USM3D 

Full N-S, finite 
volume 

Mac 64 bit 

2 CPUs 

16 hr 

1.9 M cells for half 
model without T-tail 


TRANAIR 

Full potential + 
viscosity 
(boundary 
layer) 

Linux 

Workstation 

1CPU 

2h, 28min 

1.7M cells for full 
model 
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Fig. 1 Comparison of flight measured and calculated (CFD) pressure on Hyper- X/X-43 vehicle 



Fig. 2 Grumman Gulfstream El (GEI) business jet. 
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(a) Domain Discretization 



(b) Surface Mesh 

Fig. 3 Aerodynamic model of the GUI aircraft wing. 
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Fig. 6 C p plot at span station 385. 



(a) C p distribution on wing surface 



(c) C p at station 230 


(b) C p at station 145 



(d) C p at station 385 


Fig. 7 Typical C p plots at various locations 
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